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ABSTRACT 

Gravitational instability (GI) of a dust-rich layer at the midplane of a gaseous circumstellar disk 
is one proposed mechanism to form planetesimals, the building blocks of rocky planets and gas giant 
cores. Self-gravity competes against the Kelvin- Helmholtz instability (KHI): gradients in dust content 
drive a vertical shear which risks overturning the dusty subdisk and forestalling GI. To understand 
the conditions under which the disk can resist the KHI, we perform three-dimensional simulations 
of stratified subdisks in the limit that dust particles are small and aerodynamically well coupled to 
gas, thereby screening out the streaming instability and isolating the KHI. Each subdisk is assumed 
to have a vertical density profile given by a spatially constant Richardson number Ri. We vary Ri 
and the midplane dust-to-gas ratio po and find that the critical Richardson number dividing KH- 
unstable from KH-stablc flows is not unique; rather Ri CT it grows nearly linearly with /j,q for /iq = 
0.3-10. Plausibly a linear dependence arises for /jq <C 1 because in this regime the radial Kepler 
shear replaces vertical buoyancy as the dominant stabilizing influence. Why this dependence should 
persist at /iq > 1 is a new puzzle. The bulk (height-integrated) metallicity is uniquely determined 
by Ri and fiQ. Only for disks of bulk solar metallicity is Ri CI n ~ 0.2, close to the classical value. 
Our empirical stability boundary is such that a dusty sublayer can gravitationally fragment and 
presumably spawn planetesimals if embedded within a solar metallicity gas disk ^4x more massive 
than the minimum-mass solar nebula; or a minimum-mass disk having ^3x solar metallicity; or some 
intermediate combination of these two possibilities. Gravitational instability seems possible without 
resorting to the streaming instability or to turbulent concentration of particles. 

Subject headings: hydrodynamics — instabilities — planetary systems: protoplanetary disks — planets 
and satellites: formation 
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1. INTRODUCTION 

In the most venerable scenario for forming planetesi- 
mals, dust particles in circumstellar gas disks are imag- 
ined to settle vertically into thin sublayers ("subdisks") 
sufficiently dense to undergo gravitati onal instability 
(|Safronovlll96l IGoldreich & Wardl[l97l for a review of 
this and other ways in which planetesimals may form, 
see Chiang & Youdin 2010, hereafter CY10). Along 
with this longstanding hope comes a longstanding fear 
that dust remains lofted up by turbulence. Even if 
we suppose that certain regions of the disk are de- 
void of magnetized turbulence because they are too 
poorly ionized to sustain m agnetic activity (jGammiel 
Il996t iBai k, Goodman! l2009f) . the dusty sublayer is sus- 
ceptible to a Kelvin- Helmholtz shearing instability (KHI; 
IWeidenschiliml [19801 Fl 
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1.1. Basic Estimates 

The KHI arises because dust-rich gas at the midplane 
rotates at a different speed from dust-poor gas at alti- 
tude. The background radial pressure gradient dP/dr 
causes dust-free gas at disk radius r to rotate at the 
slightly non-Keplerian rate 



where Qk is the Kepler angular frequency, 
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is a dimensionless measure of centrifugal support by pres- 
sure, and p e is the dens ity of gas (e.g., iNakagawa et al.1 
119861 : ICuzzi et"aLl 119931 ). The numerical evaluation is 
based on the minimum-mass solar nebula derived by 
CY10. Unlike dust-free gas, dust-rich gas is loaded by 
the extra inertia of solids and must rotate at more nearly 
the full Keplerian rate to remain in centrifugal balance. 
Variations in the dust-to-gas ratio pa/Pg with height z 
result in a vertical shear dv^ / dz from which free energy 
is available to overturn the dust layer. 

The shearing rate across a layer of thickness Az is given 
to order of magnitude by 
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where Pd/p g = Po at the midplane and Pd/Pg <C 1 at 
altitude (for more details see CY10, or § $TT\&2\ of this 
paper) . For po ^> 1 the velocity difference Av<p saturates 
at a speed nVL^r ~ 25(r/AU) 1 / 14 m s^ 1 , well below the 
gas sound speed c s ~ lkms -1 . That the flow is highly 
subsonic motivates what simulation methods we employ 
in our study. 

We might expect the flow to be stabilized if the Brunt- 
Vaisala frequency 
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of buoyant vertical oscillations is much larger than the 
vertical shearing rate. For the order-of-magnitude eval- 
uation in Q we approximate the vertical gravitational 
acceleration g as the vertical component of stellar grav- 
ity — il^-Az (no self- gravity), and the density gradi- 
ent p- l dpjdz ~ (p d + p g ) _1 A(p d + p g )/Az ~ (p d + 
p g ) _1 Apd/Az. The last approximation relies in part on 
the dust density pd changing over a lengthscale Az much 
shorter than the gas scale height. Both \dv^/dz\ and 

^Brunt snrm k as po decreases 

For two-dimensional, heterogeneous, unmagnetized 
flow, a necessary but not sufficient condition for insta- 
bility is given by the Richardson number: 
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< 1/4 is necessary for instability 

(5) 



(|Mileslll961t see the textbook bv lDrazin fc Reidl I2004D 
The Richardson number is simply the square of the ratio 
of the stabilizing Brunt frequency (U]) to the destabiliz- 
ing vertical shearing frequency ([3]). The critical value 
of f/4 arises formally but can al so be derived heurist i- 
cally by energy arguments (e.g., iChandrasekharl Il981h . 
The Richardson criterion does not formally apply to our 
dusty subdisk, which represents a three-dimensional flow: 
the KHI couples vertical motions to azimuthal motions, 
while the Coriolis force couples azimuthal motions to ra- 
dial motions. (For how the Ric hardson criterion ma y 
not apply to magnetized flows, see lLecoanet et aLll2010L ) 
Nevertheless we may hope the Richardson number is use- 
ful as a guide, as previou s works have assumed (iSekival 
[1991: lYoudin fc Shull200l lYoudin fe Chiang|l2004D ~ 

In this spirit let us use the Richardson criterion to es- 
timate the thickness of a marginally KH-unstable dust 



6 But not indefinitely. In the limit fio - > 0, the vertical shearing 
and Brunt frequencies reach mini ma set by pressure and t emper- 
ature gradients in gas (see, e.g., Knobloch & Spruit 1985). The 
limit fio — > is not relevant for our study and not captured by 
either Q or (gj. 
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Since the gas scale height H g = c s /T2k and r\ ~ (H g /r) 2 , 
equation ([5]) indicates that for po > 1 the marginally 
unstable dust sublayer is ^Ri x / 2 H g /r ~ 0.02i?i 1//2 times 
as thin as the gas disk in which it is immersed. Those 
KH-unstable modes that disrupt the layer should have 
azimuthal wavelengths — and by extension radial wave- 
lengths, because the Kepler shear turns azimuthal modes 
into radial ones — that are comparable to Az. Shorter 
wavelength modes cannot overturn th e layer, while longer 
wave length modes grow too slowly (|G6mez fc Ostrikerl 
l200l. 

How does disk rotation affect the development of the 
KHI? In a linear analysis. Ilshitsu fc Sekival (| 2003f) high- 
light the role played by the Keplerian shear, character- 
ized by the strain rate 
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in limiting the growth of KH-unstable modes. The radial 
shear is implicated because azimuthal motions excited by 
the KHI are converted to radial motions by the Coriolis 
force; moreover, the non-axisymmetric pattern excited 
by the KHI is wound up, i.e., stretched azimuthally by 
the radial shear. The Kepler rate \dd^/dlar\ is at least 
as large as WBrunt, and can dominate the latter when 
po is small. This suggests that Ri does not capture all 
the relevant dynamics — a concern already clear on formal 
grounds. In this paper we address this concern head-on, 
using fully 3D numerical simulations to assess the role of 
the Richardson number in governing the stability of the 
dust layer. 

1.2. Our Study in Relation to Previous Numerical 
Simulations 

Three-dimensional shearing box simulations of the KHI 
in dusty subdisks, performed in the limits that dust is 
perfectly coupled to gas and disk self-gravity is negli- 
gible, demonstrate the importance o f the Kepler shear. 
Compared to rigidly rotati ng disks (|G6mez fc Ostrikerl 
120051 : Ijohansc n et al.|[2006l). radially shearing disks are 
far more stable (|Chiang|l200l IBarrancol [20091 ). The rel- 
evance of Ri, or lack thereof, may be assessed b y simu- 
lating flows with initially spatially constant Ri ([Sekival 
199& lYoudin fc Shu! [2002). and varying Ri from run to 
run to see whether dust layers turn over. IChianel ([2008, 
hereafter C08) found that when p > 1, dust layers for 
which Ri < 0.1 overturn, while those for which Ri > 0.1 
do not. In retrospect, we might have anticipated this 
result, that the critical value Ri C nt dividing stable from 
unstable runs lies near the canonical value of 1 /4, at least 

7 This order-of-magnitude expression for the dust layer thick- 
ness, and the related equation I0 which approximates the ver- 
tical shear, are each smaller than their counterparts given by 
lYoudin & Shu (2002, page 499, first full paragr aph) by a factor 
of (1 + /1). This is because Youdin & Shu (2002) evaluate quanti- 
ties deep inside the layer, within a density cusp at the midplane, 
whereas we are interested in quantities averaged across the entire 
layer. The difference does not change cither our conclusions or 
theirs. 
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for /iQ > 1, because in this regime of parameter space all 
the frequencies of the problem are comparable to each 
other: \dv$jdz\ ~ WBrunt ~ \dQ%/dlnr\ ~ fix when 
Ho > 1 and i?z « 0.1-1. But other simulations of C08 
also make clear that Ri does not alone determine stability 
under all circumstances. For /io « 0.2-0.4, i?i cr ;t was dis- 
covere d to drop subs tantially to ^0.02 (see his runs S9- 
S12). Chiang (200$) speculated tha t the baroclinic na- 
ture of the flow may be responsible (|Knobloch fc S pruit 
I1985L I1986T ). but no details were given. 

In addition to being left unexplained, the findings 
of C08 require verification. Parameter space was too 
sparsely sampled to discern trends with confidence. Con- 
cerns about numerics — e.g., biases introduced by box 
sizes that were too small, resolutions too coarse, and runs 
terminated too early — also linger. At least one numerical 
artifact marred the simulations of C08: the KHI mani- 
fested first at the "co-rotation" radius where the mean 
azimuthal flow speed was zero (see his figure 8). But in a 
shearing box, by Galilean invariance, there should be no 
special radius. It was suspected, but not confirmed, that 
errors of interpolation associated with the grid-based ad- 
vection scheme used by C08 artificially suppressed the 
KHI away from co-rotation. 

For the problem at hand, t he spectral code developed 
bvlBarranco fc Marcus! (|2006l ) and modified bv lBarrancol 
(|2009t hereafter B09) to treat mixtures of dust and gas is 
a superior tool to the grid-based ZEUS code utilized by 
C08. Working in Fourier space rather than configuration 
space, the simulations of the KHI by B09 did not be- 
tray the co-rotation artifact mentioned above. Spectral 
methods, often used to model local (WKB) dynamics, are 
appropriate here because the structures of interest in the 
subdisk have dimensions tiny compared to the disk radius 
(by at least a factor Ri 1 / 2 !] according to equation[6]) and 
even the gas scale height. At the same computational ex- 
pense, spectral algorithms typically achieve greater effec- 
tive spatial resolution than their grid-based counterparts 
(|Barranco fc Marcus! 120061 ). Another advantage enjoyed 
by the B09 code is that it employs the anelastic approxi- 
mation, which is designed to treat subsonic flows such as 
ours. Having filtered away sound waves, anelastic codes 
are free to take timesteps set by how long it takes fluid to 
advect across a grid cell (which themselves move at the 
local orbital velocity in a shearing coordinate system). 
By contrast, codes such as ZEUS take mincing steps lim- 
ited by the time for sound waves to cross a grid cell. The 
latter constraint is the usual Courant condition for nu- 
merically solving problems in compressible fluid dynam- 
ics. It was unnecessarily applied by C08 to a practically 
incompressible flow. 

In this paper we bring all the advantages of the spec- 
tral, anelastic, shearing box code of B09 to bear on the 
problems originally addressed by C08. We assess nu- 
merically the stability of flows characterized by constant 
Richardson number Ri, systematically mapping out the 
stability boundary in the parameter space of Ri, mid- 
plane dust-to-gas ratio /io, and bulk metallicity Ed/S g 
(the height-integrated surface density ratio of dust to 
gas). Though our simulations may still be underresolved, 
we rule out box size as a major influence on our results. 
We offer some new insight into why Ri is not a sufficient 
predictor of stability. And in the restricted context of 



our constant Ri flows, we assess the conditions neces- 
sary for the midplane to become dense enough to trigger 
gravitational instability on a dynamical time. 

1.3. The Perfect Coupling Approximation vs. The 
Streaming Instability vs. Turbulent Concentration 
Between Eddies 

Following C08 and B09, we continue to work in the 
limit that dust is perfectly coupled to gas, i.e., in the 
limit that particles are small enough that their frictional 
stopping times t st0 p in gas can be neglected in compar- 
ison to the dynamical time f^ 1 - The perfect coupling 
approximation allows us to screen out the streaming in- 
stability which relies on a finite stopping time and which 
is most powerful when particles ar e marginally coupled, 
i.e., when t s = Sljcistop ~ 0.1-1 (jYoudin fc Goodman! 
2005). Numerical simulations have shown that when 
an order-unity fraction of the disk's solids is in parti- 
cles having t s = 0.1—1, the streaming instability clumps 
them strongly and paves the way for gra vitational in- 
stability (e.g., QohanieniEaD HQ03, [20091). The parti- 
cle sizes corresponding to t s = 1 depend on the prop- 
erties of the background gas disk, as well as on the 
particle's shape and internal density; under typical as- 
sumptions, marginally coupled particles are decimeter to 
meter-sized. 

It remains debatable whether a substantial fraction of 
a disk's solid mass is in marginally coupled particles at 
the time of planetesimal formation, as current propos- 
als relying on the streaming instability assume. Particle 
size and shape distributions ar e not well constrain ed by 
observations (though see, e.g., iWilner et al.l 120051 who 
showed that centimeter-wavelength fluxes from a few T 
Tauri stars are consistent with having been emitted by 
predominantly centimeter-sized particles). Measuring r s 
in disks also requires knowing the gas density, but direct 
measurements of the gas density at disk midplanes do not 
exist. Marginally coupled particles — sometimes referred 
to as "meter-sized boulders" — also face the longstanding 
problem that they drift onto the central star too quickly, 
within hundreds of ye ars from distanc e s of a few AU in 
a minimum-mass disk. Uohansen et al.l (|2007t ) claimed to 
solve this problem by agglomerating all the boulders into 
Ceres-mass planetesimals via the streaming instability 
before they drifted inward. Their simulation presumed, 
however, that all of the disk's solids began boulder-sized. 
The concern we have is that even if particle-particle stick- 
ing could grow bould ers (and sticking is exp ected to stall 
at centimeter sizes; IBlum fc WurmI 120081: CY10), the 
disk's solids may not be transformed into boulders all 
at once. Rather, marginally coupled bodies may initially 
comprise a minority population on the extreme tail of 
the particle size distribution. Unless they can transform 
themselves from a minority to a majority within the ra- 
dial drift timescale, they would be lost from the nebula 
by aerodynamic drag. 

By focussing on the dynamics of the smallest, most 
well entrained particles having r s -C 1, our work com- 
plements that which relies on the streaming instability. 
We would argue further that the well coupled limit is po- 
tentially more relevant for planet formation. If even the 
smallest particles having sizes <C cm can undergo gravita- 
tional collapse to form kilometer-sized or larger planetes- 
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imals, nature will have leapfrogged over the marginally 
coupled regime, bypassing the complications and uncer- 
tainties described above. 

Particle clumping is not restricted to marginally cou- 
pled particles via the streaming instability. Small r s 
particles also clump within the interstices of turbulent , 
high vorticity edd ies (|Maxevlll987HEaton fc Fesslerll 19941: 
ICuzzi et"aLll20Q8l and references therein; for a review, 
see CY10). This particle concentration mechanism pre- 
sumes some gas turbulence, which may be present in 
the marginally KH-unstable state to which dust settles. 
Our simulations cannot capture this phenomenon. How- 
ever, on the scales of interest to us, turbulent clumping 
might only be of minor significance. Particles of given 
istop are concentrated preferentially by eddies that turn 
over on the same timescale. Thus the degree of con- 
centration depends sensitively on particle size and the 
turbulent spectrum. At least in Kolmogorov turbulence, 
the smallest eddies concentrate particles most strongly 
because they have the greatest vorticity. The small- 
est eddies at the inner scale of Kolmogorov turbulence 

have sizes l\ ~ i/ 3 ^ 4 ^ 4 /Svl^ 2 , where v is the molecu- 
lar kinematic viscosity, and t a and 5v are the turnover 
time and velocity of the largest, outer scale eddy. Given 
6v ~ r]fl K r ~ 25(r/AU) 1 / 14 m/s, t ~ ft" 1 , and val- 
ues of v based on the nebular model of CY10, we esti- 
mate that £i ~ 10 3 (r/AU) 127 / 56 cm. This is far smaller 
than the sublayer thicknesses Az ~ 0.02i?i 1 / 2 _ff g ~ 

2 x 10 9 (ffi/0.1) 1/2 (Y/AU) 9/7 cm considered in this pa- 
per. Moreover, the lifetimes of the particle clumps on a 
given eddy length scale should roughly equal the eddy 
turnover times, which for the smallest eddies are of or- 
der U ~ y/T^/Sv ~ 10 2 (r/AU) 55 / 28 s. We do not ex- 
pect such rapid fluctuations in particle density, occur- 
ring on such small length scales, to affect significantly 
the evolution of the slower, larger scale KHI. Turbulent 
clumping may only serve as a source of noise on tiny 
scales. The possibility that turbulent clumping could 
still b e significant on larger scales is still bein g investi- 
gated (|Hogan fc Cuzzi|[2007t ICuzzi et al.|[2008h . 

The perfect coupling approximation prevents us from 
studying how particles sediment out of gas into dusty 
sublayers, but it does not stop us from identifying what 
kinds of sublayers are dynamically stable to the KHI. A 
subdisk with a given density profile is either dynamically 
stable or it is not, and we can run the B09 code for 
many dynamical times (typically 60 or more) to decide 
the answer. In a forthcoming paper we will combine the 
B09 code with a settling algorithm that will permit us to 
study how dust settles from arbitrary initial conditions, 
freeing us from the assumption that the density profile 
derives from a constant Richardson number. 

1.4. Organization of this Paper 

Our numerical methods, including our rationale for 
choosing box sizes and resolutions, are described in ^J2j 
Results are presented in £j3]and discussed in £|4] 

2. METHODS 

The equations solved by the B09 code are rederived in 
£12.11 Initial conditions for our simulations are given in 
£12.21 The code itself is briefly described in £12.31 Our 
choices for box size and resolution are explained in £12.41 



2.1. Equations 

The equations we solve are identical to equations (12a- 
e) of B09. We outline their derivation here, filling in steps 
skipped by B09, adjusting the notation, and providing 
some clarifications. This section may be skimmed on a 
first reading. 

We begin with the equations for an ideal gas perfectly 
coupled to pressureless dust in an inertial frame: 



dv 
It 



VP 



Pd + Pg 



dt 

d(Pd/p e ) 
dt 



= 0, 
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0) 
(10) 

(11) 
(12) 



where d/dt is the convective derivative, p g (d) is the den- 
sity of gas (dust), P is the gas pressure, and T is the 
gas temperature. Under the assumption that they are 
perfectly coupled, gas and dust share the same velocity 
v, and the dust-to-gas ratio is conserved in a Lagrangian 
sense. The background potential is provided by the cen- 
tral star of mass M: $ = —GM/\/r 2 + z 2 , where r is the 
cylindrical radius and z is the vertical distance above the 
disk midplane. There are five equations for the five flow 
variables v, p g , p<i, P, and T. The thermodynamic con- 
stants include the specific heat Cy = 3?/ (7 — 1) at con- 
stant volume, the ideal gas constant 5i = Cp — Cy, the 
specific heat Cp at constant pressure, and 7 = Cp/Cy. 
Equation (|11[) is equivalent to the condition that the flow 
be isentropic [d{Pp~ 1 )/dt — 0]. The code which solves 
the fluid equations actually employs an artificial hyper- 
viscosity to damp away the smallest scale perturbations 
f£ )2.3[) : in writing down equations (|5))- (fT2"j) . we have omit- 
ted the hyperviscosity terms for simplicity. 

We move to a frame co-rotating with dust-free gas at 
some fiducial radius r — R. This frame has angular fre- 
quency ilp given by © with ft K = (GM/R 3 ) 1 / 2 . We 
define a velocity v ma x using the pressure support param- 
eter 77, as given by ([2]): 



fmax = V\ r =R ^kR- 



(13) 



The velocity v max is the difference in azimuthal velocity 
between a strictly Keplerian flow and dust-free gas; it is 
the maximum possible difference in velocity, attained at 
large po, between gas at the midplane and gas at altitude. 
The quantities v max , ?7, and the background radial pres- 
sure gradient are equivalent; specifying one specifies the 
other two. Our numerical models are labeled by u max - 

In addition to moving into a rotating frame, we also 
replace the usual cylindrical coordinates (r, <fr, z) with lo- 
cal Cartesian coordinates x = r — R, y — (<p — i7pt)P, 
and z@ Keeping terms to first order in \x\ ~ \z\ ~ r/R 
(see the discussion surrounding equation[S]) and dropping 



8 Throughout this paper we alternate freely between subscripts 
(x,y,z) and (r, </>, z). 
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curvature terms, the momentum equation ([5} reads 



~dt 



-2^kzx v+3r2 K a;x— fl K zz- 



Pa 



-VP-^f^r/Px 



(14) 

where d/dt = d/dt + Vid/dxi (i — x,y,z). On the 
right-hand side, the first term is the Coriolis accelera- 
tion, the second combines centrifugal and radial grav- 
itational accelerations, the third represents the vertical 
gravitational acceleration from the star, and the last term 
arises from the centrifugal acceleration in a frame rotat- 
ing at J7f 7^ ^k- The remaining fluid equations appear 
the same as (|9|)- (fT2|) . except that v is now measured in 
a (rigidly) rotating frame. 

We measure all flow variables relative to a time- 
independent reference state (subscripted "ref"): 



v ro f 



P = P 



rcf 



= Pg,re£ + Pg 



T = T ref + T 



Pa = Pd.rct + Pa = Pa ■ 

The reference state is defined as follows. It is dust-free 
(pa,ref = 0) and has constant gas temperature T ro f. The 
gas in the reference state does not shear, either in the 
radial or vertical directions, but rotates with a fixed an- 
gular frequency Of in the inertial frame (hence v ro f = 
in the rotating frame). In the reference state there exists 
a radial pressure gradient directed outward 



' ' y/ '" !f = 2tlfaR = 2n K v B 



2;, rcf 



(15) 



and a vertical pressure gradient balanced by vertical tidal 
gravity 



1 dR 



rcf 



-rcf 



dz 



(16) 



Equation (fTBf together with equation ([12")) and the as- 
sumption of constant T ro f implies that the reference gas 
density p gjre f and pressure P re f have Gaussian vertical 
distributions in z with scale height H g = \/dtT rc f /Qk- 
For simplicity we neglect the radial density gradient 
(dpg^ e f/dr = 0), as did B09. This reference state should 
not be confused with our equilibrium states of interest 
( §2.2)1 . which do shear and which do contain dust. The 
reference state merely serves as a fiducial. 

The flows of interest are subsonic. Mach numbers e = 
v/c s peak at v max /c a ~ Cs/(QkP) ~ 0.02 for gas sound 
speeds c s ~ 1 km/s at R ~ 1 AU. Such flow is nearly 

incompressible: \p g \/p e , Ie { ~ \P\/P Ic f ~ \T\/T Ief ~ e 2 . 
Invoking the anelastic approximation, we keep only terms 
leading in e in any given equation. Equations (J9)), (fTU|) . 



and (fT2"j) reduce to 
dp 



dt 



PgV ' V = lit + V ' (PgV) W V ' (( ° s ' rcf V) = ° 

(17) 

d(pa/Pg) _ rf(pd/p g ,rcf) _ dp = Q 
(it c?t 



P 



= /i = 



Pg,rof 



Pg.ref 



-3?T rof + 3fiT (19) 



where we define p = pa/ Pg.xcf = Pa/Pg,rtf and the 
pressure- like enthalpy h = P/p g , rc f, and henceforth for 
convenience drop all tildes on pa, p, and v (but not the 
other variables related to gas). The rightmost equal- 
ities of ([T7]). ((15)1 . and match equations (12b), 
(12c), and (12e) of B09. The anelastic approxima- 
tion has b een employed in the study of atmo spheric 
convection (IQgura fc Phillips! 119621: IGoughl(T969l). stars 
(|Gilman fc Glatzmaierlll98lD. and vortices in proto plan- 
etarv disks (|Barranco fc Marcus! [2000l [20051 l2006h . By 
eliminating the time derivative in the continuity equation 
(fT7|) . we effectively "sound-proof" the fluid. The simula- 
tion timestep is not limited by the sound-crossing time 
but rather by the longer advection time. 

We rewrite our energy equation f|l 1 1) as follows: replace 
-V • v with d\\\p g /dt = -dhxT/dt + dlnP/dt to find 
that 



C P 



dT 1 dP 



dt 



dt 



-V ■ VP ref 



-rcf 



i-v • (2^k?7Px + 0|zz) 



(20) 



where for the second line we dropped dP/dt in com- 
parison to v • VP re f, and for the third line we replaced 
p~] e{ VP re f using (TT51) and (fTO|) . Equation (f2"0|) matches 
(12d) of B09 except that for the right-hand side he has 

a coefficient equal to 1 + T/T rc {, which we have set to 
unity. 

Finally, to recover the form of the momentum equation 
(12a) of B09, first consider the pressure acceleration and 
isolate the contribution from dust-free gas (— p K T 1 VP): 



1 



-VP = - 



Pa 



Now expand 
1 



- — ) VP - —VP 



Pa 



Pe Pi 

Ivp^ 



-VP. 



(21) 



-VP 



Py 
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where for the last line we used (fT5|) and (IT51) . Insertion 
of ([21]) and ([22]) into (| 14[) yields the anelastic momentum 
equation (12a) of B09: 

dv T 

— = -2fi K z x v + 3QkXX + — — ^fi^iix + f^zz) - V/i 
at T rc f 



fi 



fx + 1 



1 + — ) (2^?yi?x + Q^zz) - Vh 
T lc f / 



(23) 



which isolates the driving term due to dust. 



2.2. Initial Conditions 
Equilibrium initial conditions (superscripted "f") are 
specified by five functions: p = $ , T = T^, ft. = /i^, 
Pg = /5g^, and v = v^. For p^ , we use flows characterized 
by a globally constant Richardson number (|Sekivalll998t 
lYoudin fc Sh"v3[200l IChiandl200ll . The conditions Ri = 
constant, dp g /dz <C dpd/dz, and g — —Q^z (no self- 
gravity) yield 



p\z) 



1 



1/2 



1. 



.l/(l + /i ) 2 + (z/z d ) 2 
where po is the initial midplane dust-to-gas ratio and 

Ri 1 / 2 v 



(24) 



fir- 



(25) 



is a characteristic dust height. The dust density peaks 
at the midplane and decreases to zero at 



Jp (2 + p ) 

z = ±z max = ± — z d (26) 

1 + Mo 

which is consistent with our order-of-magnitude expres- 
sion ©. Neither equation (|2"4"|) nor the code accounts for 
self-gravity and therefore we are restricted to modeling 
flows whose densities are less than that required for the 
Toomre parameter of the subdisk to equal unity (CY10; 
see also Sj4j. For the minimum- mass disk of CY10, this 
restriction is equivalent to p <j 30. Input model param- 
eters include po, Ri, and v max . 
For the gas, we assume 



rt = 



(27) 



(initially isothermal) and solve vertical hydrostatic equi- 
librium for (the ^-component of equation I23[) : 



dz 



(28) 



The functional form for h'(z) is not especially revealing 
and so we do not write it out here. For simplicity we 
assume that does not depend on x. From b) and 
ft = it follows from (H]) that 



Pg,refft t 



ref 



(29) 



The fractional deviations p g V/°g,ref and P* 1 /P le f from the 
reference state are very small, of order p^ {v max / c s ) 2 Ri. 



It remains to specify v^. Using the conditions on fv 
stated above, we solve for the equilibrium (steady-state) 
solution to equation ([23]) : 



vt=vt = 



tf(z) + l 



(30) 



In our reference frame rotating with the velocity of dust- 
free gas at R, the first term on the right side of ([30]) 
accounts for the standard Kepler shear, while the second 
term describes how dust, which adds to inertia but not 
pressure, speeds up the gas. 
To fv we add random perturbations 

Ap(x,y,z) = A(x,y)p^ (z)[cos(irz /2z d ) + s\\\{kz /2z d )] . 

(31) 

The amplitude A{x,y) is constructed in Fourier space 
so that each Fourier mode has a random phase and 
an amplitude inversely proportional to the horizontal 
wavenumber: A oc fcj 1 = (fc 2 . + fc 2 ) -1 / 2 . Because our 
box sizes are scaled to z max , our Fourier noise ampli- 
tudes are largest on scales comparable to the dust layer 
thickness. Thus those modes which are most likely to 
overturn the layer are given the greatest initial power. 
The perturbations are also chosen to be antisymmetric 
about the x-axis so that no extra energy is injected into 
the system. We take the root-mean-squared amplitude 
^ rms = {A 2 ) 1 / 2 of the perturbations to be 10~ 4 or 10~ 3 . 

In summary, three input parameters po, Ri, and v max 
determine our isothermal equilibrium initial conditions 
(equations [Ml HH and |3U[) O The equilibrium solution 
for p(z) is then perturbed (equation l31[) by a root-mean- 
squared fractional amount ^4 rms . The parameters of pri- 
mary interest are po and Ri. For the remaining parame- 
ters v max and A rms we consider three possible combina- 
tions: (v max , A rms ) = (0.025c s , 10~ 4 ) for our standard 
runs; (0.025c s , 10~ 3 ) to probe larger initial perturba- 
tions; and (0.05c s , 10~ 4 ) to assess the effect of a stronger 
radial pressure gradient. 

Note that specifying po and Ri (and t> max , though this 
last variable is fixed for all of our standard runs) specifies 
the entire dust and gas vertical profiles, Pd{z) and p g (z), 
and by extension the bulk height-integrated metallicity, 
Sd/Sg = J pddz/ J p g dz. We do not give an explicit 
expression for Sd/S g because it is cumbersome and not 
particularly revealing. The bulk metallicity is in some 
sense the most natural independent variable because its 
value is given by the background disk (for ways in which 
the bulk metallicity may change, e.g., by radial particle 
drifts, see CY10). We will plot our results in the space 
of /io, Ri, and Ed/^g, keeping in mind that only two of 
these three variables are independent. 

2.3. Code 

We use t he spectral, anelastic , shea ring box code de- 
veloped by iBarranco fc Marcusl (|2006[ ) and modified by 

9 While our initial conditions are isothermal, the temperature of 
the flow can change because of adiabatic compression/expansion 
and because our artificial hyperviscosity dissipates the highest 
wavenumber disturbances. These temperature changes are frac- 
tionally tiny because the flow is highly subsonic. 
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B09 to simulate well-coupled gas and dust. The code em- 
ploys shearing periodic boundary conditions in r, peri- 
odic boundary conditions in <f>, and closed lid boundaries 
in z; the vertical velocity v z is required to vanish at the 
top and bottom of the box (z = ±L z /2). 

Spectral methods approximate the solution to the fluid 
equations as a linear combination of basis functions. The 
basis functions describe how the flow varies in space, and 
the coefficients of the functions are determined at every 
timcstep. For each of the periodic dimensions, a stan- 
dard Fourier basis is used, while for the vertical direc- 
tion, Chebyshev polynomials are employed. Whereas in 
r and grid points are evenly spaced, the use of Cheby- 
shev polynomials in z has the effect that vertical grid 
points are unevenly spaced; points are concentrated to- 
wards the top and bottom boundaries of the box, away 
from the midplane where the dust layer resides. Thus to 
resolve the dust layer vertically, we need to increase the 
number of vertical grid points N z by an amount dispro- 
portionately large compared to the numbers of radial and 
azimuthal grid points N r and N$. See Ej2.4l for further 
discussion. 

Spectral codes have no inherent grid dissipation; en- 
ergy is allowed to cascade down to the smallest resolved 
length scales through nonlinear interactions. To avoid 
an energy "pile-up" at the highest wavenumbers, we dis- 
sipate e nergy using an artific i al hyp erviscosity, given in 
§3.3.3 of lBarranco fc Mareusl (pOOl ). 

Simulations satisfy the Courant-Fricdrichs-Lewy 
(CFL) condition which states that the CFL number, 
defined as the code timestep divided by the shortest 
advection time across a grid cell, be small. In the 
shearing coordinates in which the code works, that 
advection time is the cell dimension divided by the local 
velocity over and above the Keplerian shear, i.e., orbital 
velocities are subtracted off before evaluating the CFL 
number. All simulations reported in this paper are 
characterized by CFL numbers less than about 0.1. 

2.4. Box Size and Numerical Resolution 

Our standard box dimensions are (L r ,L^,L z ) = 
(6.4, 12.8, 8)z max and the corresponding numbers of grid 
points are (N r ,N,p,N z ) = (32,64,128). By scaling our 
box lengths Li to z max and fixing the numbers of grid 
points Ni, we ensure that each standard simulation en- 
joys the same resolution (measured in grid points per 
physical length) regardless of Ri, /i U5 and u max - The ver- 
tical extent of the dust layer between z = ±z max is re- 
solved by 22 grid points (this is less than [128/(8z max )] x 
2^max = 32 because the Chebyshev-based vertical grid 
only sparsely samples the midplane). The radial and 
azimuthal directions are resolved by 10 grid points per 
%z max length. We choose our resolution in the vertical 
direction to be greater than that of the horizontal direc- 
tions because the dust layer has finer scale structure in 
z: the dust layer becomes increasingly cuspy at the mid- 
plane as no increases. We prescribe the same resolution 
in the radial and azimuthal directions (L^/N^ = L r /N r ); 
experiments with different resolutions in r and <p gener- 
ated spurious results. 

Too small a box size can artificially affect the stability 
of the dust layer, because a given box can only support 
modes having integer numbers of wavelengths inside it. 
Small boxes may be missing modes that in reality over- 
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Fig. 1. — Testing box sizes at fixed numerical resolution. 
For our standard box, (L r ,L^,L z ) = (6.4, 12.8, 8)z m ax and 
(N r ,N,p, N z ) = (32, 64, 128). In each panel we vary one box dimen- 
sion while keeping the other two dimensions fixed at their standard 
values. In the top panel we vary L z at fixed resolution N z /L z . In 
the middle and bottom panels, and L r are varied in turn. All 
simulations in this figure have fio = 10, Ri = 0.1, f ma x = 0.025e s , 
A rms = 10~ 4 , and use code units p g!lc f(z = 0) = fix = Hg = 1- 
Doubling the box dimensions from our standard values changes 
when the average vertical kinetic energy peaks by only a few orbits 
at most. The average () is performed over all r and at fixed 
2 = 0. 

turn the layer. We verify that for all runs in which the 
dust layer overturns, the KH mode that most visibly dis- 
rupts the layer spans more than one azimuthal wave- 
length. Typically 3-5 wavelengths are discerned across 
the box. 

To more thoroughly test our standard choices for 
Li, we study how systematic variations in box length 
affect how the instability develops. For this test, 
we adopt a fixed set of physical input parameters, 
(Ri, Ho, Uma x) = (0.1, 10, 0.025c s ), which should lead to 
instability (jChiand 120081) . Our diagnostic is the time 
evolution of the vertical kinetic energy at the midplane: 
{n{t)v1 (t))/2, where the average is over all r and 4> at 
fixed z — and time t. We vary Li and Ni in tan- 
dem to maintain the same resolution from run to run, 
thereby isolating the effect of box size. Figure Q] shows 
how doubling one of the box dimensions while fixing the 
other two alters the time history of (/zi^)/2. Panel (a) 
demonstrates that our standard choice for L z = 8z max 
is sufficiently large because the curves for L z — 8,z max 
and L z = 16z max practically overlap. Panels (b) and (c) 
show that our standard choices for L^ = 12.8z max and 
L r = 16z max are somewhat less adequate. The peak of 
the curve for (L^, Nq) — (12.Sz max , 64) is delayed by two 
orbits compared to that for (L^,N^,) = (25.6z max , 128), 
and the curve for (L r , N r ) — (6Az max , 32) peaks an orbit 
earlier than that for (L r ,N r ) = (12.8z max , 64). Never- 
theless these time differences are small compared to the 
total time to instability, about 10 orbits. Moreover, the 
errors point in opposite directions. Thus we expect our 
choices for L^ and L r to partially compensate for each 
other so that any error due to our box size in calculating 
the time to instability will be less than ~1 orbit. 



We test how robust our results are to numerical resolu- 
tion by re-running a few simulations at twice the normal 
resolution (doubling N{ while fixing L,). Results at high 
resolution are given in ^3-31 Every simulation is run for 
at least ten orbits. A typical run performed at our stan- 
dard resolution takes approximately 2.5 wall-clock hours 
using 56 processors on the Purdue Steele cluster. A high- 
resolution run takes about 32 wall-clock hours. 

3. RESULTS 

In our standard simulations, we fix w max and A lms 
while systematically varying Ri and fj-o from run to run. 
Our systematic variations of Ri and /io correspond to 
systematic variations in £d/£ g ; recall that only two of 
the three parameters Ri, /io, and £d/S g are indepen- 
dent. For each /io 6 {0.3, 1, 3, 10} we adjust Ri until the 
threshold value Ri C mt dividing stable from unstable runs 
is determined to within 0.1 dex. 

Deciding by numerical simulation whether a given dust 
layer is stable or not is unavoidably subject to the finite 
duration of the simulation. We define our criteria for 
deciding stability in ^3.11 Results are given in ^3.21 and 
tested for robustness in £13.31 

3.1. Criteria for Stability 

Stability is assessed by two quantities: the midplane 
vertical kinetic energy 

(nv 2 z )/2 as a function of t 

where the average is performed over r and <f> at fixed 
z = and t, and the dust density profile 

(fj,) as a function of z and t 

where the average is performed over r and <fi at fixed 
z and t. By definition, in an "unstable" run, (fiv 2 )/2 
grows exponentially over several orbital periods, and (/i) 
deviates from its initial value (fr) by more than 15%. 
"Stable" simulations satisfy neither criterion. Some runs 
are "marginally unstable" in that they satisfy the first 
but not the second criterion. At the end of the standard 
ten-orbit duration of a marginally unstable run, we find 
the kinetic energy continues to rise, suggesting that were 
the run to be extended for longer than ten orbits, the 
dust layer would eventually overturn. In every instance 
where we extend the duration of a marginally unstable 
run, we verify that this is the case. Thus "marginally 
unstable" is practically synonymous with "unstable." 

Examples of unstable and stable runs are shown in Fig- 
ure [2j In the unstable simulation, after t«6 orbits, the 
kinetic energy rises exponentially. At t w 9 orbits, the 
dust layer overturns and the midplane dust-to-gas ratio 
falls by more than 60%. By contrast, in the stable sim- 
ulation, after an initial adjustment period lasting ^3 or- 
bits during which the midplane value of (/x) decreases by 
10%, the kinetic energy drops by orders of magnitude to 
a nearly constant value and shows no evidence of further 
growth. 

Figure [3] shows the evolution of (A — r, 0, z) and 

(n(z)) for the same unstable run of Figure^] The velocity 
data are sampled at a single (x, y) position at the center 
of our simulation box. The radial and vertical velocities 
| IV | and \v z \, initially zero, grow to become comparable 



with the shearing velocity \v t j > \. Figure U displays corre- 
sponding snapshots of (j,(y, z), taken at a single radius x 
near the center of our box. Though the data in Figures 
[3] and |4] are sampled at particular radial locations in our 
box, we verify that the instability develops similarly at 
all locatio ns — as it should — unlike the ZEUS-based sim- 
ulations of lChiand (|200l . 

3.2. Stability as a Function of Ri, /io, and £d/£ g 

Figure [5] maps the stable and unstable regions in 
(Ri, no) space, for fixed v max = 0.025c s and A rms = 10~ 4 . 
Figures[6]and[7]portray the same data using alternate but 
equivalent projections of parameter space: (Ri, Sa/^g) 
and (/io,Sd/Sg), respectively. 

These plots demonstrate that there is no unique value 
of fficrit- Rather i?i cr it is a function of /io, or equivalently 
a function of Ed/S g . For bulk metallicities Ed/£ g near 
the solar value, Ri cr it is found to be close to the classical 
value of 1/4. But as Sd/E g decreases below the solar 
value, fficrit shrinks to ~0.01 or even lower. A least- 
squares fit to the four midpoints (evaluated in log space) 
in Figure [5] dividing neighboring stable points (in black) 
and unstable points (in red or red outlined with black) 
yields i?i cr it oc /ij'°- This same fit, projected into metal- 
licity space, is shown in Figures [6] and [7j in metallicity 
space the stability boundary is not a power law. 

As Figure [7] attests, dust-to-gas ratios /io as high as 
~8 can be attained in disks of solar metallicity without 
triggering a shear instability: see the intersection be- 
tween the dashed curve fitted to our standard resolution 
data, and the dotted line representing solar metallicity. 
This intersection occurs at /io ~ 7. Were we to re-fit 
the dashed curve using the higher resolution data repre- 
sented by triangles, the intersection with solar metallicity 
would occur at /io closer to 8. 

A dust-to-gas ratio of /io ~ 8 is within a factor of 
~4 of the Toomre threshold for gravitational fragmen- 
tation in a minimum-mass disk (CY10; @. We can 
achieve the Toomre threshold by simply allowing for a 
gas disk that is ^4x more massive than the minimum- 
mass nebula. Alternatively we can enrich the disk in 
metals to increase Ed/S g . Extrapolating the boundary 
of stability (dashed curve) in Figure [7] to higher Ed/S g 
suggests that the Toomre threshold /io ~ 30 could be 
achieved for minimum-mass disks having ^3x the solar 
metallicity. The sensitivity to metallicity is also exem- 
plified by Figure [21 For the same /io = 10, the dust layer 
based on a near-solar metallicity of Ed/S g = 0.013 over- 
turns, whereas one derived from a supersolar metallicity 
of Sd/S g = 0.030 remains stable. 

3.3. Tests at Higher Resolution, Higher j4 rms , and 

Higher u max 

We test how robust our determination of Ri cr it is to 
numerical resolution by redoing our simulations for fj-o — 
0.3 and 10 with double the number of grid points in each 
dimension. The results are overlaid as blue triangles in 
FiguresEl El and[7] At /io = 0.3, increasing the resolution 
does not change Ri C rit from its value of 0.009. At /io = 
10, Riant shifts downward from 0.3 to 0.2. Although 
we have not strictly demonstrated convergence of our 
results with resolution, and although high resolution data 
at other values of /io are missing, it seems safe to conclude 
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Fig. 2. — Sample unstable (top) and stable (bottom) dust lay- 
ers. In the unstable case, the layer overturns and mixes dust-rich 
gas with dust-poor gas, causing the dust-to-gas ratio at the mid- 
plane to drop by a factor of ~3 after 10 orbits (top left). As the 
instability unfolds, the vertical kinetic energy amplifies exponen- 
tially from t « 5-10 orbits (top right). At fixed the layer is 
stabilized by increasing the Richardson number or equivalently the 
height-integrated metallicity S^/Sg. In the stable run, the dust 
profile changes by less than 15% (bottom left) while the kinetic 
energy, after dropping precipitously, shows no indication of grow- 
ing (bottom right). The two runs shown use v max = 0.025c a and 
A rms — 10 

that the slope of the stability boundary in Ri-^lq space 
is close to, but decidedly shallower than, linear. 

We also test the sensitivity of our results to A Yma . In- 
creasing A rms by an order of magnitude to 10 ~ 3 shifts 
-R*crit upward by <; 0.2 dex at /io < 1, but leaves i?i cr it 
unchanged at larger /i (Figure [5]). B09 also reported 
some sensitivity to A lms . 

Tests where u max was doubled to 0.05c s reveal no 
change in Ri cl it (data not shown). 

4. SUMMARY AND DISCUSSION 

Where a protoplanetary disk is devoid of turbulence in- 
trinsic to gas, dust particles settle toward the midplane, 
accumulating in a sublayer so thin and so dense that 
the dust-gas mixture becomes unstable. If the first in- 
stability to manifest is self-gravitational, dust particles 
are drawn further together, possibly spawning planetes- 
imals. If instead the layer is first rendered unstable by 
a Kelvin- Helmholtz-type shearing instability (KHI), the 
resultant turbulence prevents dust from settling further, 
pre-empting gravitational collapse. In this paper we in- 
vestigated the conditions which trigger the KHI, hoping 
to find a region of parameter space where the KHI might 
be held at bay so that planetesimals can form by self- 
gravity. 

A fundamental assumption underlying our work is 
that turbulence intrinsic to gas can, in some regions of 
the disk, be neglected. There is some consensus that 
near disk midplanes, in a zone extending from ~1 to 
at least ^10 AUs from the parent star, gas may be 
too poorl y ionized to susta in ma gnetohydrodynam ic tur- 
bulence Aligner fc Nelson] [2006t iBai fc Goodman! [200l 
iTurner et al.l 12010). Presumably if the magnetorota- 
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Fig. 3. — Snapshots of absolute values of the three velocity 
components (top panels) and horizontally averaged dust-to-gas 
ratio (bottom panels), both as functions of height, at three in- 
stants in time. For this unstable run, (Ri, fio, Umax, A rms ) = 
(0.1, 10, 0.025e s , 10~ 4 ). Velocities are taken from a grid point near 
the middle of the box. The vertical shear dv^/dz inside the dust 
layer weakens with time as dust is more uniformly mixed with gas, 
and as the radial and vertical velocities grow at the expense of the 
azimuthal velocity. 



tional instability (e.g.. iBalbusI 1200 9) cannot operate at 
the midplane, disk gas there is laminar — pending the un- 
certain ability of magnetic ally active surface layers to 
stir the disk interior (e.g., ITurner et al.ll2010[ ). or the 
discovery of a p urely hydrodynamic form of turbulence 
(Lithwick 2009). To get a sense of how laminar disk 
gas must be to permit dust sublayers to form, Chiang & 
Youdin (2010) compared the height to which dust parti- 
cles are stirred in an "alpha" -turbulent disk to the thick- 
ness of the sublayer ([6|). They estimated that the for- 
mer is smaller than the latter when the dimensionless 
turbulent diffusivity a^3x 10- 4 fl K t stop (r/AU) 4/7 for 

£stop < fijf . To place this requirement in context, a 
values for magnetically active zones are typically quoted 
to be greater than ^10~ 3 . Whether magnetically dead 
zones are sufficiently passive for dust to settle into sub- 
layers remains an outstanding question. 

Modulo this concern, we studied the stability of dust 
layers characterized by spatially constant Richardson 
numbers Ri using a th ree-dimensional, spectral, anelas- 
tic, shearing box code (|Barranco fc Marcusl[20CMl 
that models gas and dust as two perfectly coupled fluids 
(|Barrancoll2009[ ). We found that stability is not charac- 
terized by a single critical Richardson number. Rather 
the value of i?i C rit distinguishing layers that overturn 
from those that do not is a nearly linear function of 
the midplane dust-to-gas ratio /xo (Figure E]) ■ Dust-rich 
sublayers having /^o ~ 10 have i?? C rit ~ 0.2 — near the 
canonical value of 1/4 — while dust-poor sublayers hav- 
ing fiQ ss 0.3 (still orders of magnitude dustier than well- 
mixed gas and dust at solar abundances) have Ri cl n as 
low as 0.009. 

Pre vious studies (e.g.. iSekival Il998t lYoudin fc Shu! 
I^OllYoudin fc Chiandl20o¥ assumed a universal criti- 
cal Richardson number of 1/4. This popular assumption 
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Fig. 4. — Snapshots of fi(y,z), sampled at r = R (x = 0; the 
central slice of the simulation box) for the same unstable run 
shown in Figure [3] The box size parameters are (L r , La, L z ) = 
(0.05, 0.1, 0.063)-ff g , larger than what is shown in the figure, which 
zooms in for more detail. 

seems correct only for dust-rich layers having /io so large 
they are on the verge of gravitational instability. For less 
dusty midplanes, the assumption appears to be incorrect. 
O ur numerical r esults are roughly consistent with those 
of [Chiang (200§), who also found evidence that i2i cr it de- 
creases with decreasing /io. Comparing his Table 2 with 
our Figure [5] shows that his constraints on Ri CT it are, for 
the most part, compatible with those presented here, for 
the range /io ~ 0.3-10 where ou r respective da ta overlap. 
Our findings supersede those of lChiana (|2008T ) insofar as 
we have explored parameter space more finely and sys- 
tematically, at greater and more uniform resolution, with 
numerical methods better suited for subsonic flows. 

Our results turn out to be consistent with the classical 
Richardson criterion — which states only that Ri < 1/4 is 
necessary, not sufficien t, for i n stabil ity — even though the 
criterion as derived bv I Miles! (| 1 9 6 lh applies only to two- 
dimensional flows, which our dust layers are not. Our 
simulations demonstrate that the criterion can still serve 
as a useful guide for assessing stability in disks having 
bulk metallicities ranging from subsolar to slightly super- 
solar values — with the proviso that the actual Richard- 
son number dividing KH-stable from KH-unstable flows, 
while < 1/4, is generally not equal to 1/4. 

Why isn't the Richardson criterion for instability suf- 
ficient in rotating dust disks? The criterion consid- 
ers the competition between the destabilizing vertical 
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Fig. 5. — Mapping the boundary of stability in the space of initial 
Ri and no. Red points correspond to unstable dust layers, whose 
dust-to-gas ratios (fi) change by more than 15%, and whose verti- 
cal kinetic energies grow exponentially, within the 10-orbit duration 
of the simulation. Black points mark stable dust layers satisfying 
neither criterion. Red points outlined in black signify marginally 
unstable layers, whose kinetic energies rise but whose dust-to-gas 
ratios change by less than 15%; these are essentially equivalent to 
red points without outlines, because every marginally unstable run 
that we extend beyond 10 orbits eventually becomes fully unstable. 
Runs performed at twice the standard resolution appear as trian- 
gles. Downward pointing triangles symbolize stable runs, upward 
triangles are unstable, and upward pointing triangles in black out- 
line are marginally unstable. All simulations U.SG Arms — 10 -4 and 
» ma i = 0.025c s . There is no unique value for the critical Richard- 
son number separating stable from unstable dust layers. Rather, 
a least-squares fit to the data from our standard resolution runs 
yields Ri CI a oc /i 10 , shown as a dashed line. The classical bound- 
ary Ri CIl t = 0.25 is plotted as a dotted line. 

shear and the stabilizing influence of buoyancy, which 
causes fluid parcels to oscillate about their equilibrium 
positions at the Brunt- Vaisala frequency. However, 
there exists another stabilizing influence, ignored by the 
Richardson number, pr ovided by the radial Kepler shear 
(jlshitsu fc Sekiva1 [2003'). In the limit /io <C 1, the Brunt 
frequency (@| becomes negligible relative to the Kepler 
shearing frequency (J7J , suggesting stability now depends 
on the competition between the destabilizing vertical 
shear and stabilizing radial Kepler shear. We expect 
the flow to be stable as long as the Kepler shear can 
wind up unstable eigenmodes to higher radial wavenum- 
bers before their amplitudes grow large enough to trig- 
ger nonlinear effects. This suggests that we replace the 
Richardson number with a "shearing number," defined 
by analogy as the square of the ratio of the Kepler shear- 
ing frequency to the vertical shearing frequency: 



Sh: 



\dn/dlnr\ 

(dv^/dz) 2 



oc 



Az 

Av<p 



oc Ri- 



1 



Mo 



Mo 



(32) 



where we have used (|3|) and ([6]). By assuming Sh is 
constant for marginally stable dust profiles, we arrive at 
the relation 



fficrit oc Mo for Mo < 1 



(33) 



What is surprising is that this trend, although expected 
to hold only for /io -C 1, appears to hold approximately 
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Fig. 6. — Mapping the boundary of stability in the space of initial 
Ri and bulk (height-integrated) dust-to-gas ratio Ed /Eg. The data 
are identical to those in Figure [5] The labeling convention is also 
the same, except that the triangles representing high-resolution 
runs have adjusted their orientation so that they point towards 
the stability boundary. The same least-squares fit from Figure[5]is 



projected here as a dashed curve. Solar metallicity E^/Eg 



: 0.015 



Pro. 

(Loddcrs 2003) is indicated by a dotted line. The critical value 
Ri CT li dividing stable from unstable dusty subdisks trends with 
metallicity. This trend was only hinted at in the data of C08. 
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Fig. 7. — Mapping the boundary of stability in the space of mid- 
plane dust-to-gas ratio /io and bulk (height-integrated) dust-to-gas 
ratio E^/Eg. The data are identical to those in Figure [5] The la- 
beling convention is also the same, except that the triangles repre- 
senting high-resolution runs have adjusted their orientation so that 
they point towards the stability boundary. The same least-squares 
fit from Figure[5]is pr ojected here as a dashed curve. Solar metal- 
licity E(j/Eg = 0.015 (Loddcrs 2003) is indicated by a dotted line. 
A minimum-mass solar nebula requires ~ 30 for gravitational 
instability to ensue on a dynamical time (CY10). Extrapolating 
the boundary of stability to fj,o ss 30 suggests that metallicitics 
roughly ^3 times solar would be required for dynamical gravita- 
tional instability in a minimum-mass disk. The required degree of 
metal enrichment would be proportionately less in more massive 
disks. 
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Fig. 8. — How the stability boundary changes with stronger initial 
perturbations. This figure is the same as Figure [5] except that all 
data correspond to A rms = 10 -3 . For comparison with /Irms — 
10" 4 , the same best-fit line of Figure [5] is reproduced here. Not 
much changes, except that ffi cr it shifts upward by 0.2 dex at /iq = 
0.3. 

for all /io, according to our simulation results in Figure [5j 
For ytto ;> 1, we would have expected from (|32|) that ffi cr it 
asymptote to a constant; but it does not. Our higher res- 
olution runs do suggest the stability curve slightly flat- 
tens at fj,o « 10, but such deviations seem too small to 
be fully explained using arguments relying purely on the 
shearing number. 

To explai n the observed t r end, we might co-opt the 
methods of llshitsu fc Sekival (|2003D . who linearized and 
numerically integrated the 3D equations of motion for 
the dust layer. For their particular choice of background 
vertical density profile, they solved for the maximum 
gro wth factors for the most unstable KH modes (see 
also lKnobloch fc SpruitllT985l who considered the axisym- 
metric problem). We would need to replace their as- 
sumed profile with our profiles having spatially constant 
Ri. Perhaps our numerically determined stability curve 
i?« cr it(Sd/Sg) corresponds to a locus of fixed maximum 
growth factor. 

Gravitational instability occurs on a dynamical time 
when the dust layer's Toomre Q ~ M/\2irr 3 p s (l + Mo)l 
reach es unity (jToomrd Il964t IGoldreich fc Lvnden-Belll 
1965). For p g given by the minimum- mass solar neb- 
ul a, this occurs when pp ~ 30, fairly independently of 
r (jChiang fc Youdinll2010[ ). Of course in more massive 
gas disks (greater p g ), the requirement on po is pro- 
portionately lower. Figure [7] shows that for disks hav- 
ing bulk metallicities Sd/S g equal to the solar value of 
0.015, the dusty sublayer can achieve po ~ 8 before 
it becomes KH unstable. Taken at face value, such a 
marginally KH-stable subdisk, embedded in a gas disk 
having 30/8 ss 4 times the mass of the minimum- mass 
solar nebula, would undergo gravitational instability on 
the fastest timescale imaginable, the dynamical time. 
The case that planets form from disks several times more 
mass ive than the minimum-mass solar nebula is p lausible 
(e.g., IGoldreich etaLll2004t iLissauer et aLll2009f) . 

An alternate way of crossing the Toomre threshold is 
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to allow the bulk metallicity E^/Sg to increase above the 
solar value of 0.015. Extrapolating the boundary of sta- 
bility in Figure [7] to /io ~ 30 suggests that metallicities 
roughly ^3 times solar would be required for dynamical 
gravitational instability in a minimum-mass disk. There 
are several proposed ways to achieve supersolar metal- 
licities i n some portions of t he disk, among them radial 
pileups (lYoudin k, Shu 2002 ) or dissipative gravitational 
instabil ity (lWardlll976t]Coradini et al]H98lTlWardll200tt 
lYoudinll2005at lYoudirJ l2005bt see also the introduction 
of lGoodman fc Pindorll2000( ). 

None of the ways we have outlined for achieving grav- 
itational instability rely on the streaming instability or 
turbulent concentration of particles, mechanisms that we 
have criticized in Nevertheless our scenarios may 

be too optimistic because all our dust profiles are predi- 



cated on the assumption of a spatially constant Ri. This 
assumption tends to generate strong density cusps at the 
midplane that might not be present in reality. In a forth- 
coming paper we will relax the assumption of spatially 
constant Ri and measure the maximum fi attainable, as 
a function of metalllicity Sd/S g , by simulating explicitly 
the settling of dust towards the midplane. 
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Science Foundation, in part through TeraGrid resources 
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